Quantum Monte Carlo studies of edge magnetism in chiral graphene nanoribbons 
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We investigate chiral graphene nanoribbons using projective quantum Monte Carlo simulations 
within the local Hubbard model description and study the effects of electron-electron interactions 
on the electronic and magnetic properties at the ribbons' edges. Static and dynamical properties 
are analyzed for nanoribbons of varying width and edge chirality, and compared to a self-consistent 
Hartee-Fock mean-field approximation. Our results show that for chiral ribbons of sufficient width, 
the spin correlations exhibit exceedingly long correlation lengths, even between zigzag segments 
that are well separated by periodic armchair regions. Characteristic enhancements in the magnetic 
correlations for distinct ribbon widths and chiralities are associated with energy gaps in the tight- 
binding limit of such ribbons. We identify specific signatures in the local density of states and low- 
energy modes in the local spectral function which directly relate to enhanced electronic correlations 
along graphene nanoribbons. These signatures in the local density of states might be accessed by 
scanning tunneling spectroscopy on graphene nanoribbons. 
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I. INTRODUCTION 

Graphene nanoribbons (GNRs), laterally confined, 
nanometer wide one-dimensional long carbon strips, are 
currently intensively examined with respect to their elec- 
tronic properties and potential for future graphene-based 
electronic devicesP^ Various fabrication strategies have 
been explored recently in order to achieve high-quality 
GNRs, such as the unzipping of carbon nanotubegSHS, 
or the direct chemical synthesis of GNRsP^ Transport 
measurements report sizable band gaps in GNRsplwhich 
in general depend on the edge geometry. In case of GNRs 
with armchair edges (aGNRs), these gaps are a conse- 
quence of the lateral quantum confinement, while for 
the highly symmetric GNRs with zigzag edges (zGNRs) , 
an energy gap is theoretically predicted to arise due to 
the spontan eous sp in polarization established along the 
zigzag edges IHiSllfil Such edge magnetism stems from the 
electron-electron interactions in the high density of elec- 
tronic states at the Fermi level which is predominantly 
localized near the zigzag edgesP^ Across the ribbon, the 
two edges are magnetically anti-aligned to each otheil^, 
thus respecting Lieb's theorem of zero net magnetiza- 
tion for these bipartite lattice structures. ^ Ab-ini tio 
density functional theory (DFT) calculation ^ 1 ' 16 ' 19 ' 20 1 , as 
well as mean-field theory (MFT) based on the Hartree- 
Fock decoupling within a tight-binding Hubbard model 
description of zigzag GNRs predict long-ranged ferro- 
magnetic order along the edges of zigzag GNR in the 
ground stateP^ Exceedingly large magnetic correlation 
lengths were indeed observed in unbiased numerical stud- 
ies of the edge magnetism using quantum Monte Carlo 
(QMC) simulations^, the density matrix renormaliza- 
tion group^l and exact diagonalizationPSHU Steps to- 
wards a rigorous proof of the emergence of magnetic mo- 
ments on zigzag edges have been put forward in Ref. 1261 




FIG. 1. (Color online) (a) Armchair graphene nanoribbon 
of width w = 12, denoted 12-aGNR. (b) Chiral ribbon of 
chirality (3,1) (8 = 13.9°), and width w = 6, denoted as 
6(3, l)-cGNR. (c) Zigzag nanoribbon of width w = 6, denoted 
6-zGNR. In all cases, the two sublattices of the bipartite 
lattice structure are shown using light and dark spheres and 
the ribbons' unit cells are indicated by dashed lines. Armchair 
(zigzag) segments along the upper edge within the unit cell 
in all cases are highlighted in red (blue). 



In the more generic case of chiral GNRs (cGNRs), charac- 
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terized by the presence of both zigzag and armchair seg- 
ments, the presence of low-energy edge states for edges 
with sufficiently long zigzag contributions has been in- 
vestigated previouslyp2H22l Such ribbon geometries are in 
particular relevant for recent experiments that provided 
direct evidence for the presence of GNR edge statesP^ES] 

Recent theoretical investigations of chiral GNRs based 
on Hubbard-model MFT and ab-initio DFT studies pro- 
mote the existence of spin polarized edges, similar to 
the scenario in zigzag GNRp^SI. Here, we apply un- 
biased large-scale QMC simulations to study both the 
static and dynamical properties of chiral GNRs to assess 
the robustness of the edge magnetism with respect to 
the armchair segments, and compare our QMC results to 
self-consistent MFT calculations. 

The remainder of this paper is organized as follows: 
The next section introduces the nomenclature of the dif- 
ferent GNR geometries that we analyze. After a short 
overview of the numerical method that we employ, in 
Sec. Ill we first focus on the static magnetic proper- 
ties and discuss the spin correlations for GNRs of dif- 
ferent chiralities. In Sec. IV we study the local density 
of states (LDOS), which we calculate from the single- 
particle Green's function. We identify low-energy fea- 
tures in the spectral functions and local variations in the 
LDOS characteristic for correlated edges, which may be 
probed via scanning tunneling spectroscopy. 



II. RIBBON GEOMETRIES AND METHODS 

In the following, we consider several different GNR ge- 
ometries, and use the standard notations to specify their 
edge structure: For armchair GNRs, we denote a ribbon 
with w dimer lines parallel to the ribbon direction as w- 
aGNR, e.g., Fig.[]|a) shows a 12-aGNR. A zigzag GNR 
ribbon with w zigzag lines is denoted as w-zGNR. The 
6-zGNR structure is illustrated in Fig. [ljc). For chiral 
GNRs, wc follow the notation in Ref. (311 wherein the 
direction of the edge is specified by a translation vector 
(n, m). In particular, the edge of a (n, m) ribbon consists 
of a repeated unit of m armchair units followed byre-m 
zigzag units. The width of the chiral GNRs in terms of 
the number w of parallel zigzag units between the edges is 
indicated by the notation w(n,m)-cGNR. The example 
in Fig. |TJb) shows a 6(3, l)-cGNR. While the armchair 
case corresponds to (n,m) — (1,1) (more specifically, a 
w(l, l)-cGNR is identical to a 2u>-aGNR) and the zigzag 
case to (n,m) — (1,0), we employ the more convenient 
notation introduced above to highlight these symmetric 
edges. The chirality of a GNR edge can also be charac- 
terized by the chirality angle 9, i.e., the angle between 
the translation vector and the closest zigzag line, which 
is obtained from the identity 

sine = J 3 - (V^ (1) 

V 4 \n 2 + nm + m 2 J w 



In chiral ribbons it spans the range 0° < 6 < 30°, where 
for zigzag edges 8 — 0°, and for armchair edges 6 = 30°. 
The distance between nearest neighbor sites is denoted 
here by clq and in the following serves as the unit for all 
spatial distances. In particular, the length of the trans- 
lation vector (n,m), i.e. the extent of the ribbon's unit 
cell, is given by a = aoVn 2 + nm + m 2 . In all cases, we 
employed periodic boundary conditions (PBC) along the 
extension of the ribbons. 

The influence of electron-electron interactions and the 
resulting magnetic properties of chiral GNRs in its most 
basic form is captured by the Hubbard model 

H = - t Y. fcW + hx ") + U E n *T "4 • (2) 

Here, t denotes the hopping between nearest-neighbors 
and U the onsite repulsion. The operator cJ CT cre- 
ates a spin-er fcrmion on site i and the local density is 
defined by n ia = c\ a c ia . Estimates of the local Coulomb 
repulsion in graphene- based m aterials put the ratio U/t 
in the region near unity^EIES i n the following, we thus 
consider the weak coupling regime U/t < 2, where bulk 
magnetic order is absenpSHS] Most of the simulation 
results are shown for U/t = 2, where the magnetic prop- 
erties are more pronounced and thus may be robustly 
detected within finite size QMC simulations. Further- 
more, we consider here the case of half-filling, where the 
total number of electrons equals the number of lattice 
sites. In this case, sign-problem free QMC simulations 
of the Hamiltonian H can be performed on all the above 
introduced topologies. 

We use a projective determinantal QMC method, 
which allows to extract ground state expectation values 
of an arbitrary observable O (we set h = 1) 

(tto|Q|tto) (y T \e-° H Oe-° H \V T ) 
<*o|*o> e>^> (* T |e- 2eH |* T ) ' 

by projection from a trial wave function I^t), which 
is taken here to be the eigenstates of the free system 
(U = 0). The projection parameter is chosen suffi- 
ciently large, such that convergence to the ground state 
|*o) is guaranteed. Depending on the detailed GNR 
structure, values of 6 between 6 = 60/t and 100/i are 
required to ensure convergence. We use a symmetric 
Suzuki- Trotter decomposition with an imaginary time 
discretization of At = 0.05/t, such that discretization 
errors are well below the size of the statistical errors. 
Furthermore, a SU(2) symmetric Hubbard-Stratonovich 
decoupling of the Hubbard interaction was employed in 
order to ensure the explicit conservation of spin rota- 
tional symmetry. Details of the projective QMC method 
can be found in Ref. |33J We simulate lattices of Nq unit 
cells along the ribbon and impose PBC. For all the con- 
sidered ribbons, we specify the value of Nq available for 
the largest simulated structure in detail below. As an ex- 
ample, for the 6(3, l)-cGNR structure shown in Fig.[ljb), 
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FIG. 2. (Color online) Real-space spin correlations illustrated 
by disks proportional to (Si • Sj) across GNRs of various chi- 
ralities of width w = 6 (corresponding to a 12-aGNR in the 
armchair case) at U/t — 2. In each case, the reference site 
is indicated by an arrow. Red (blue) circles correspond to 
positive (negative) values. 



which contains 48 sites per unit cell, we simulated peri- 
odic ribbon rings with up to Nc = 12 unit cells, i.e. a 
total of 576 sites. 

In addition to the QMC simulations we solve 
the self-consistent mean-field equations for the 
magnetization, which we obtain from the Hartree- 
Fock decoupling of the Hubbard interaction 
Hit "4 n it (njj.) + (nif)nn - (n it )(n 4 ), ignoring 
the fluctuation term. For regular ribbon geometries, 
a Fourier transformation along the ribbon direction is 
performed. 



III. MAGNETIC CORRELATIONS 

As a consequence of the finite density of states at the 
Fermi level, cGNRs have been predicted to exhibit an 
instability toward a spin-polarized edge upon the intro- 
ductions of electron-electron interactions. Based on MFT 
this edge magnetism is expected to persist essentially 
over the full rang e of c hiralities, but for the near armchair 
limit of 9 = 30° I n the zigzag limit, for which the 
edge polarization shows a maximum, exceedingly large 
spin correlation lengths have been confirmed by means 
of unbiased QMC simulations.^ In contrast, GNRs with 
perfect armchair edges of widths w ^ 3k + 2 (with in- 



teger k = 0,1,2,..) yield exponentially decaying (short 
range) spin correlations. This can be directly related to 
the semiconducting nature of such GNRs. Armchair rib- 
bons with w = 3fc + 2 in the non interacting limit are 
metallic and exhibit unique behavior, discussed in the 
following. 

Due to the presence of repeated armchair segments in 
chiral GNR, the issue arises, whether the long-range spin 
order that results within the MFT approximation is ac- 
tually robust to quantum fluctuations. To address this 
issue, we performed QMC simulations for GNRs of dif- 
ferent chiralities and widths. To probe for tendencies 
towards edge magnetism within our QMC simulations, 
we measure the equal-time spin-spin correlations between 
sites and in particular monitor the spin correlation func- 
tion C(r) = (So • S r ) as a function of the distance r along 
the edges of the GNRs. 

The general correlation pattern along the ribbon edge, 
as well as between the two edges and into the bulk, is 
illustrated in the real-space representations of the corre- 
lations in Fig. [2]for three characteristic cases, also con- 
sidered in Fig. [3] The reference sites are indicated by 
arrows. The extremely short range correlations in the 
armchair GNR in Fig. [2la) steadily increase from the 
6(3,l)-cGNR in Fig. gbjto the zigzag GNR in Fig.^c). 
In all cases the correlations decay quickly into the bulk 
and are anti-aligned to each other on opposite ribbon 
edges, respecting the sublattice structure. Like in zigzag 
GNRs, the correlations in the chiral GNR are largest 
along the zigzag segments interrupted by armchair el- 
ements where correlations are suppressed. Compared 
to the zigzag case, the correlations into the bulk decay 
quicker and are strongest along the direction towards the 
closest opposite edge. 

In Fig. [3j we show results for the correlation function 
C(r) as a function of the distance r from a given refer- 
ence site for GNRs of different chiralities and width w = 6 
(in the armchair case, this corresponds to a 12-aGNR in 
standard notation). These results were obtained on the 
longest ribbons of each chirality type accessible in our 
QMC simulations. The correlations are shown among 
edge sites of the same sublattice with respect to a refer- 
ence site centrally located within the zigzag segment for 
chiral GNRs. Consecutive data points correspond in each 
case to a spatial separation of a. For better comparison 
among the different cases, we plot all distances r in units 
of ao, the nearest neighbor distance on the graphene lat- 
tice. The inset of Fig. [3] illustrates this procedure for the 
case of a 3,l)-cGNR. Note that for armchair and zigzag 
GNRs, all outer edge sites provide equivalent reference 
points. 

For this intermediate ribbon width (w — 6), we observe 
the clear tendency towards quasi long-range spin corre- 
lations. Starting from the nearly exponentially decaying 
correlations in the aGNR, they become increasingly dom- 
inated by power law behavior with increasing length of 
the zigzag segments, before essentially long-range spin 
correlations prevail. We estimate the spin correlation 
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FIG. 3. (Color online) Spin correlation function C(r) along 
the edge of GNR of various chiralities of width w — 6 (corre- 
sponding to a 12-aGNR in the armchair case). For each case, 
the employed number of unit cells JVc is indicated. Lines 
show the fitting function Cfit(r) for each GNR case. The in- 
set in the lower left part of the figure indicates the position of 
the reference site for the case of a (3,l)-cGNR edge (arrow), 
as well as the distance to the corresponding sites (triangles) 
within further unit cells, each separated by a distance a. 
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FIG. 4. (Color online) Spin correlation function C(r) along 
the edge of GNR of various chiralities of width w = 12 (corre- 
sponding to a 24-aGNR in the armchair case). For each case, 
the employed number of unit cells JVc is indicated. Lines show 
the fitting function Cfu(r) for each GNR case. 



length £ for the various chiralities by fitting a function 
that combines both, exponential and power-law behavior: 



Cla t (r) oc r-"e- r /« + (L — r )-*e-( L - r M 



(4) 



Here the second term accounts for the PBC, with L de- 
noting the linear extent of the ribbon along its edgel^S 
For these ribbons, the fitted values of the exponent rj 
range between 0.2 and 1.9 and we concentrate on dis- 
cussing the correlation length £ in the following. For 
the 12-aGNR armchair ribbon we extract an expectedly 
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FIG. 5. (Color online) Spin correlation function C(r) along 
the edge of aGNRs of different widths uu. Lines in this figure 
are guides to the eye. 



short correlation length of £ = 3.2(2)ao- The 6(2,1)- 
cGNR, 6(3,l)-cGNR and 6(4,l)-cGNR show a £/a which 
steadily increases from 6.5(5) via 13.6(6) to 26(2), re- 
spectively (numbers in brackets quantify the uncertainty 
in the last digit). The resulting fitting functions Cfi t (r) 
are shown along with the QMC data in Fig. [3| In the 
6(3,2)-cGNR, with two armchair segments separating the 
zigzag regions, we find a stronger exponential decay with 
£ = 5.7(l)<Zo- Similar to the armchair case this ribbon 
lacks spin polarized edges even in MFT below U/t = 1.94 
(not shown). In the zigzag ribbon, the spin correlation 
length extracted from the fit exceeds the system size even 
on this largest available system, which is in accordance 
with previous QMC results.^ 

The generally rapid decay of the spin correlations along 
armchair edges is presented in Fig.[5j which shows results 
for armchair GNRs of different widths. A noticeable ex- 
ception from this trend is provided by the 8-aGNR: In 
contrast to other armchair GNRs, where £ < 6ao, we ex- 
tract much stronger correlations. This behavior traces 
back to the metallic nature of the 8-aGNR in the U = 
tight binding limit, which is among the series w = 3fc + 2 
(with integer k = 0,1,2,...). It is not related to the 
emergence of edge magnetism, in particular since the low- 
energy states in this case are not edge-localized. Instead, 
the behavior here is similar to the Hubbard model on 
the one-dimensional chain (to which indeed the armchair 
GNR degrades for w = 2), for which U > triggers an 
instability towards a quasi-long-ranged-ordered bulk an- 
tiferromagnetic spin density wave state. 

In increasingly wider chiral GNRs, the tendency to- 
wards edge magnetism becomes more established. This 
can be seen by comparing the results in Fig. [3] for w = 6 
to those in Fig. |1J which provides data on different rib- 
bons twice the width (w = 12 for cGNR, corresponding 
to w = 24 for the armchair case). While the w = 24 
armchair ribbon still shows a fast decay of the corre- 
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FIG. 7. (Color online) Map of the single particle excitation 
gap for the U — tight-binding limit for w(n, l)-cGNR of 
varying n (i.e. chirality) and width w. 



IV. DYNAMIC SPECTRAL FUNCTIONS 

The equal-time spin-spin correlations, examined in 
Sec. Ill, allow to quantify the strength of the magnetic 
correlations along the GNRs. We next consider the dy- 
namical single-particle spectral function, which is partic- 
ularly important, as it may be directly compared with 
scanning tunneling spectroscopy experiments. The spin 
resolved local and momentum-resolved spectral functions 
are given at frequency u) by 



A itT (u>) = K*o|cfcr|*»>| 2 S(uj + E - E n ) 



(5) 



FIG. 6. (Color online) Spin correlation function C(r) along 
the edge of (a) (2,l)-cGNR and (b) (3,l)-cGNR of different 
widths. Lines in this figure are guides to the eye. 



and 



A aa (q,uj) =J2 l(*o|ca^|*„)| 2 8(u + E -E n ) , (6) 



lations, resulting in a similar correlation length as for 
w = 12, we find significantly enhanced correlations for 
the 12(2,l)-cGNR, and the 12(3,l)-cGNR. Their esti- 
mates for the £ exceed the available system sizes. En- 
hanced correlations for specific widths can also be ob- 
served in chiral ribbons. For example, they are present 
along the 4(2,l)-cGNR and the 4(3,l)-cGNR (cf. Fig. [6}. 
Compared to neighboring values of w, such enhanced cor- 
relations at w = 4 may be associated with a vanishing 
(for the 4(2,l)-cGNR) or very low (for the 4(3,l)-cGNR) 
single particle energy gap in the tight binding limit. Sim- 
ilar behavior we thus expect to also be observed for other 
specific widths of chiral ribbons, and is reflected by the 
horizontal stripes for w — 4, 7, 10, 13 in the single parti- 
cle gap map for U = shown in Fig. [7J These ribbons 
are special examples of the metallic or almost metallic 
(MAM) points identified in Ref. 1271 



respectively, where c aqa = J2j=i e I9aj Cj( QJ -) CT , with 
i(a,j) denoting the site index of the lattice site located 
at position a within the j-th unit cell, and n enumer- 
ates the full set of eigenstates \^ n ) with energy E n of 
the MFT Hamiltonian. Because of particle-hole sym- 
metry, we need to consider only the range of positive 
uj > 0. Within MFT, broken spin rotational symme- 
try leads to distinct spectral functions for both spin 
orientations, which we average to Ai(ui) — h\Ai^(ui) + 
Aji(w)]. In the QMC simulations, SU(2) symmetry is 
preserved and therefore Ai(uj) = A^{uj) — Aj^(uS), The 
Dirac <5-functions in the above formula were subjected 
to a Lorentzian broadening of Au; = 0.02t within the 
MFT calculations. Within the QMC method we can- 
not directly access real time correlation functions, and 
thus measure the momentum-resolved Green's function 
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G a (q,r) = §E„<* 

time r, instead.^ The spectral function on the real fre- 
quency axis uj may be obtained by the inversion of 



Ga(q,r) 



duj e T " A a (q,u) , 



(7) 
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FIG. 8. (Color online) Momentum-resolved single particle spectral function for a 12-aGNR for edge sites (top) and bulk sites 
(bottom), for [7 = (left), and from MFT (middle) and QMC simulations (right) for U/t = 2. 
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FIG. 9. (Color online) Momentum-resolved single particle spectral function for a 6-zGNR for edge sites (top) and bulk sites 
(bottom), for [7 = (left), and from MFT (middle) and QMC Simulations (right) for U/t = 2. 
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by means of the stochastic analytic continuation 
methodP^ The local spectral function Ai(u>) is then 
calculated from the corresponding momentum-resolved 
spectral function A a (q,u>) by integrating over the mo- 
mentum q along the ribbon direction. This procedure 
ensures that features in the local spectral function re- 
lated to different momenta q are preserved by the ana- 
lytic continuation. In the following, we consider the spec- 
tral functions for sites along the ribbon edge, denoted by 
A e d ge {q, fjj) and ^4 e d ge (w), and within the center (bulk) of 
the ribbon denoted by A^^q, uj), and Ab u ik(w), respec- 
tively. For chiral ribbons the edge site refers to a site in 
the center of a zigzag segment, unless noted otherwise. 

In order to discuss the effects of electron interac- 
tions on the LDOS, it is instructive to first examine the 
momentum-resolved local spectral functions for GNRs 
of different chiralities. For this purpose, we show in 
Figs. [8lpd1 QMC results for A e d ge (q,uj) and Abulkfew), 
along with tight-binding (U = 0) as well as with MFT 
results. 

Starting from the spectrum of the I2-aGNR (cf. 
Fig. M, the U = tight-binding and the U/t = 2 MFT 
results are in fact identical. This is due of the absence of 
a finite spin polarization in the mean-field solution at this 
coupling strength. Upon comparing to the QMC spec- 
tral functions, one has to account for inherent resolution 
limitations imposed by the analytic continuation of the 
imaginary-time QMC data. While low energy features 
are usually well defined, the analytic continuation fails 
to resolve or broadens individual features at high ener- 
gies. Here, we indeed find that the low energy region of 
the QMC spectral function compares well to the MFT re- 
sults, revealing the same low energy characteristics. Both 
exhibit a similar excitation gap of A/t « 0.14. We hence 
find no significant difference to the non-interacting case. 

In contrast, in the zigzag GNR results shown in Fig. |9j 
we observe distinct spectral features related to ed ge mag- 
netism at finite U, as already reported previousl y 1 l 46 l 47 l. 
At U = 0, the spectral function traces the low energy 
edge states with a flat dispersion within the momen- 
tum range 2/3 < qa/n < 4/3. Both MFT and QMC re- 
sults show the increase of the gap, related to the onset 
of the edge magnetism and the characteristic bending 
of the low-energy edge band with a maximum intensity 
at q = ir/a. This pronounced low-energy band, located 
along the ribbon edge, results in a distinctive low en- 
ergy peak at w max « 0.2i in the local spectral func- 
tion A c( j KC (u), i.e. the LDOS, shown for the 6-zGNR 



in Fig. 12 while being absent in Abuik(w). As noted 



previously, MFT overestimates th e gap e nergy scales at 
U/t = 2 by about a factor of twd 21 * 46 ^ , in accordance 
with the results in Fig. [12] The integrated weight of the 
QMC low energy peak in A e( j go (o;) is robustly reproduced 
by MFT, with only about ten percent deviations among 



ting of the low energy band into three separate peaks 
seen in the QMC data due to finite size effects. 

After having discussed the cases of symmetric GNR 
edges, we next consider the spectral functions for chiral 
GNRs. Figure [10] shows the results for the 6(2,l)-cGNR. 
In the tight-binding limit, one identifies in A c dge(q, w) 
a dispersing low-energy band of edge states, which for 
wider ribbons eventually transforms into a flat band near 
q = it /a, thus providing a finite DOS of edge states 
at the Fermi energy. These low energy states of the 
6(2,l)-cGNR extend further into the bulk compared to 
the pure zigzag case, which relates to the incomplete 
suppression of the spectral weight near q = n/a in 
Ahuikiq, w) (compare to the zigzag case). Finite inter- 
actions lead to a sizable single particle gap within MFT, 
with a maximum of the spectral weight centered around 
q = it J a in A e( ^ ge (q,uj). Like for the zigzag case, this re- 
sults in a prominent low-energy peak in the LDOS at the 
ribbon edge, as seen for A c a gc (uj) in Fig. 
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the two approaches. The lower left panel of Fig. 12 shows 
MFT results for a finite zigzag GNR with iV c = 36 unit 
cells, i.e. for the same finite system as employed in the 
QMC simulations. This calculation reproduces the split- 



The spec- 
tral function obtained from QMC simulations, shown in 
Fig. |10| similarly exhibits a low energy mode predomi- 
nately confined to the GNR edge, however with a signif- 
icantly smaller gap. In the LDOS, this leads to the pro- 
nounced low-energy spectral weight in A ec [g e (^uj^) shown 
in Fig. [13] As for the zigzag case, the QMC result for 
the spectral function exhibits a discrete set of isolated 
peaks at low energies; these finite size effects are again 
reproduced within MFT by considering the same finite 
size system as in the QMC simulations (cf. the lower left 
panel in Fig. 13 1. While MFT overestimates the peak po- 
sition w max , like in the zigzag case, the integrated spectral 
weight within the low-energy regime in j4 cc i gc (w) never- 
theless agrees up to about a ten percent deviation among 
the two methods. 

For the 6(3,l)-cGNR and the 12(3,l)-cGNR, results for 
the momentum-resolved spectral functions on the edges 
are shown in Fig. 11 Like in the 6(2,l)-cGNR case, one 
identifies a low-energy band of dispersing edge states, 
now symmetric around q — 0. The different gap position 
results from the additional zigzag unit for this chirality 
and the accompanied band folding.^ Within MFT, a siz- 
able single particle gap results from the Hubbard interac- 
tion, with maximal spectral weight centered at q = 0. In 
the QMC data, a similar behavior is observed, albeit with 
a smaller gap. For the 6(3,l)-cGNR, we show the local 
spectral function Ai{uS) along the ribbon edge in Fig. 14 
In addition to the emergence of the low energy peak at 
w max ~ 0.08t, this representation also exhibits the modu- 
lation of the peak intensity along the majority sublattice 
sites. The low-energy peak is more prominently observed 
at the central sites of the zigzag regions and reduces in in- 
tensity approaching the armchair regions. Furthermore, 
the MFT data for both widths in Fig. [TT] display min- 
ima at non-zero q in the low-energy band of A cl \ gc {q 1 uj) 1 
which leads to a bending of the low-energy dispersion 
near q = 0. We observe such a shift in the position of the 
minimum gap away from q = and the related bending 
in the QMC data only for the 12(3,l)-cGNR. The gap 




FIG. 10. (Color online) Momentum-resolved single particle spectral function for a 6(2,l)-cGNR for edge sites (top) and bulk 
sites (bottom), for [7 = (left), and from MFT (middle) and QMC simulations (right) for U/t = 2. 




FIG. 11. (Color online) Momentum-resolved single particle spectral function along the edge sites for a 6(3,l)-cGNR (top) and 
a 12(3,l)-cGNR (bottom), for U = (left), and from MFT (middle) and QMC simulations (right) for U/t = 2. 
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FIG. 12. (Color online) Local single particle spectral function 
for a 6-zGNR for edge sites and bulk sites. The inset focuses 
in on the QMC data at low energies. 
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FIG. 13. (Color online) Local single particle spectral function 
for a 6(2,l)-cGNR for edge sites and bulk sites. The inset 
focuses in on the QMC data at low energies. 
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FIG. 14. (Color online) Local single particle spectral func- 
tion at low energies for a 6(3,l)-cGNR from QMC simula- 
tions along the ribbon edge (shown on the bottom), within 
the majority sublattice, as indicated by the arrows. 
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FIG. 15. (Color online) Position of the low-energy peak in 
the edge site LDOS of chiral graphene nanoribbons of width 
w = 6 as a function of the chirality angle 8 within MFT 
(open symbols) and from QMC simulations (filled symbols), 
both for U = t and U = 2t. The crosses indicate the single 
particle gap in the tight-binding (U = 0) limit. 



minimum in the 6(3,l)-cGNR still remains at q = 0, at 
least within the available resolution. We will comment 
on this observation further below. 

The results of the LDOS for the different w = 6 GNRs 
are summarized in Fig. [15] which shows the position of 
the low-energy peak position w max in A e d ge (w) for GNRs 
of different chiralities. Beyond the QMC and MFT re- 
sults for the ribbons discussed above, we included in 
Fig. [15] MFT results for ribbons of various other chirali- 
ties. In addition to the results at U/t — 2, we also show 
results for U/t = 1, as well as the single particle gap of 
the GNRs in the tight-binding limit, {7 = 0. An overall 
trend towards a reduction of u} max with increasing 9 is 
observed both in the MFT and QMC data points. MFT 
however, systematically overestimates u; max - an effect 
that apparently is enhanced for increasing chirality. One 
furthermore observes a sizable increase of the U = sin- 



gle particle gap due to the finite width of the w = 6 GNRs 
for chirality angles 9 beyond about 10°. For larger 9, this 
leads to a suppression of the edge magnetism already on 
the MFT level, as mentioned for the (3,2)-cGNR case in 
the previous section. For example, the (2, l)-cGNR and 
the (3, 2)-cGNR exhibit no edge magnetism within MFT 
at U/t = 1, while for U/t — 2 the edge magnetism is sta- 
ble for GNRs below 9 w 25°. This suppression leads to 
the apparent merging of the data for the U = gap and 
those for w max at, e.g., 9 sa 25° within MFT for U/t = 2. 

In Fig. [16] we show a similar plot for w = 12 chiral 
GNRs. When comparing to Fig. [15] we find that the po- 
sition of w max does not depend strongly on the ribbon 
width w; this holds at least within the chirality range, 
where the U = tight-binding gap is small. Such weak 
w-dependence was observed already previously for the 
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FIG. 16. (Color online) Position of the low-energy peak in 
the edge site LDOS of chiral graphene nanoribbons of width 
w = 12 as a function of the chirality angle 8 within MFT 
(open symbols) and from QMC simulations (filled symbols), 
both for U — t and U = 2t. The crosses indicate the single 
particle gap in the tight-binding (U = 0) limit. 



zigzag case in Ref. [5TJ and is also seen for the (3,1)- 
cGNRs in Fig. 11 The main difference between the two 
cases, w = 6 and w = 12, arises in the large 0-region, 
where the U = single particle gap becomes of the or- 
der w max , and hence the effective antiferromagnetic inter- 
edge coupling suppresses the formation of the edge mag- 
netism already within MFT. The intermediate ^-region 
requires a more careful analysis. Consider for example 
the case of the (3,l)-cGNR: while for w — 6, the tight- 
binding gap is already sizable, of the order of the mean- 
field gap, it is strongly suppressed for the w = 12 ribbon. 
This relates to the observation (pointed out above), that 
for the w = 12 the gap minimum in A ec j ge (9,w) shifts 
away from q = 0, leading to a bending of the low-energy 
dispersion near q = 0. For the (2,l)-cGNR, we do not 
observe a corresponding bending (at q = n/a) in the 
QMC data for A e ^ Ee {q,LS) at w = 12 (not shown), while 
in MFT such a behavior is already present on the w = 6 
ribbon (cf. Fig. [ll]). In fact, for the (2,l)-cGNR a finite 
tight-binding gap is still clearly resolved on the scale of 



shift in the gap minimum position and the bending of the 
dispersion could eventually be detected. 



V. CONCLUSIONS 

We examined the interaction-induced edge magnetism 
in chiral graphene nanoribbons based on a Hubbard- 
model description using unbiased quantum Monte Carlo 
simulations. Our results for the equal-time spin-spin cor- 
relations confirm that ribbons beyond the armchair limit 
exhibit substantial ferromagnetic correlations among the 
zigzag segments along the ribbon edges. Increasing the 
ribbon width, the effect of the antiferromagnetic inter- 
edge coupling is strongly suppressed, leading to corre- 
lation lengths that are compatible to long-range ferro- 
magnetic edge magnetism. We computed the local spec- 
tral functions related to scanning tunneling spectroscopy 
experiments, and identified a characteristic low-energy 
peak along the ribbon edge, related to the formation of 
enhanced electronic correlations. The position in energy 
of this peak is consistent with an essentially linear depen- 
dence on the interaction strength and the chirality angle, 
and shows no significant width dependence. In future 
studies, the resilience of these experimentally detectable 
features upon the introduction of realistic edge disorder 
will be investigated. Whether these features prevail as 
indication for edge magnetism beyond the mean-field ap- 
proximation, will be most efficiently studied within an ef- 
fective low-energy spin-only description of the magnetic 
correlations, which allows for the treatment of signifi- 
cantly larger systems. 
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